Mass Cytometry reveals unique phenotypic patterns associated with subclonal diversity and outcomes in multiple myeloma

Multiple myeloma (MM) remains an incurable plasma cell (PC) malignancy. Although it is known that MM tumor cells display extensive intratumoral genetic heterogeneity, an integrated map of the tumor proteomic landscape has not been comprehensively evaluated. We evaluated 49 primary tumor samples from newly diagnosed or relapsed/refractory MM patients by mass cytometry (CyTOF) using 34 antibody targets to characterize the integrated landscape of single-cell cell surface and intracellular signaling proteins. We identified 13 phenotypic meta-clusters across all samples. The abundance of each phenotypic meta-cluster was compared to patient age, sex, treatment response, tumor genetic abnormalities and overall survival. Relative abundance of several of these phenotypic meta-clusters were associated with disease subtypes and clinical behavior. Increased abundance of phenotypic meta-cluster 1, characterized by elevated CD45 and reduced BCL-2 expression, was significantly associated with a favorable treatment response and improved overall survival independent of tumor genetic abnormalities or patient demographic variables. We validated this association using an unrelated gene expression dataset. This study represents the first, large-scale, single-cell protein atlas of primary MM tumors and demonstrates that subclonal protein profiling may be an important determinant of clinical behavior and outcome.


INTRODUCTION
Multiple myeloma (MM) is the second most common hematopoietic malignancy [1,2], characterized by a clonal expansion of malignant antibody-producing post-germinal-center plasma cells (PCs) within the bone marrow (BM). Although the overall survival of patients with MM has improved over time owing to the use of novel chemotherapeutic agents, targeted therapies, and autologous stem cell transplantation [3], MM remains an incurable malignancy in most patients [4,5].
One methodology to assess single-cell proteomic heterogeneity includes mass cytometry (CyTOF) [27], which allows multiplex phenotyping of single cells providing the potential for novel protein signature identification that may be associated with chemotherapeutic resistance and/or disease outcome. CyTOF utilizes transition element isotype-tagged antibodies using time-of-flight mass spectrometry, allowing the detection of up to 40 to 50 simultaneous epitopes in single cells without significant spectral overlap. Previous studies using traditional flow cytometry have identified significant changes in the protein expression of numerous markers associated with MM disease survival, relapse, and/or drug resistance using this approach [28][29][30]. Here, we used CyTOF to analyze 34 epitopes in 49 primary MM patient samples at the single-cell level. We provide an integrated, comprehensive, single-cell atlas of cell surface molecules, transcription factors, and phosphoprotein targets that were previously individually associated with drug response or resistance in MM. We detect subclonal protein profiles that are shared among patients and associated with important clinical variables, including survival.

MATERIALS AND METHODS Patient samples
All samples were referred to the Mayo Clinic Genomics Laboratory for routine clinical testing. Forty-nine cryopreserved BM samples from patients diagnosed with MM were selected for analysis by CyTOF based primarily on the abundance of PCs in the BM and having at least 500 live light-chainrestricted PCs. This retrospective study was approved by the institutional review board of the Mayo Clinic. A minimal risk waiver for consent was obtained for this study.
Clinical data FISH, S-phase calculations, and stratification for myeloma and risk-adapted therapy (mSMART, msmart.org) assessments were performed as part of routine clinical practice. FISH analysis of immunoglobulin (cIg)-stained positive PCs studies were performed as previously described [31] using probes described in Supplementary Materials and Methods. Plasma cell S-phase data were obtained from flow cytometry, as described in [32]. Each patient's treatment response at 90 days was based on the International Myeloma Working Group uniform response [33] and further described in Supplementary Materials and Methods.

CyTOF antibody panel
Thirty-seven antibody targets (Supplementary Table 1) directed against cell surface and intracellular proteins were initially designed using the webbased panel designer software Maxpar Panel Designer (Fluidigm/Standard BioTools, South San Francisco, CA) for optimal signals, minimum background due to oxidation, isotopic purity, and sufficient sensitivity for each targeted marker. Prelabelled antibodies were purchased from Fluidigm. Purified antibodies requiring conjugation were purchased from BioLegend (San Diego, CA) or Sigma Aldrich (St. Louis, MO) and labelled using the X8 polymer Maxpar antibody conjugation kit (Fluidigm) according to manufacturer's instructions. Each antibody was validated and titrated against positive and negative control cell lines.

Sample preparation for flow cytometry and CyTOF
Total BM cells were processed using ammonium, chloride, and potassium (ACK) lysis buffer (Thermo Fisher Scientific, Waltham, MA) to lyse red blood cells and centrifuged for 5 min at 1350 RPM. The pellet was washed with PBS, centrifuged for 5 min at 1350 RPM and resuspended in Chang BMC medium containing 20% fetal bovine serum (FBS) (Irvine Scientific, Santa Ana, CA). A portion of the cells was prepared for FISH as part of routine clinical testing. The remaining cells were centrifuged at 1200 rpm for 8 min and the cell pellets were resuspended in 1.5 ml of 10% DMSO (Sigma Aldrich) in Chang BMC medium, added to a cryotube, and frozen in liquid nitrogen.
Prior to CyTOF staining, three randomly selected BM samples from the liquid nitrogen biobank were thawed and stained by flow cytometry to evaluate the abundance of total viable PCs. Each cryovial was quickly thawed and the cells were immediately washed and centrifuged at 500 x g for 5 min in 20% FBS (Thermo Fisher Scientific, Waltham, MA, USA) in RPMI1640 and 1:10,000 benzonase nuclease (Sigma Aldrich). The cell pellets were resuspended in residual RPMI solution, washed again with the same medium and resuspended in cell staining buffer including 1x PBS, 3% BSA, and 0.1% sodium azide, followed by surface staining using the following antibodies: CD38 HIT2 (Cat# 303504), CD229 HLy-9.1.25 (Cat# 326108), CD138 MI15 (Cat# 356524), and CD45 HI30 (Cat# 304026) (BioLegend). Samples were processed using a BD Accuri C6 Plus Personal Flow Cytometer (BD Biosciences, San Jose, CA) ( Supplementary Fig. 1A).
For CyTOF cell staining (reagents described in Supplementary Table 2), frozen BM cells were rapidly thawed and suspended in warmed RPMI1640 with 20% FBS containing 1:10,000 benzonase nuclease. Additional CyTOF staining details are described in Supplementary Materials and Methods. Briefly, 1-3 million cells were washed, incubated with Cell-ID and resuspended in 1x Maxpar Fix and Perm Buffer. Cells were subjected to the cell surface antibody panel, washed with Maxpar Cell Staining Buffer (CSB) and stored overnight at −80°C. Cells were subjected to the intracellular antibody panel, fixed with a fresh 1.6% formaldehyde solution (Thermo Fisher Scientific), washed with CSB and resuspended in 1 ml of Cell-ID Intercalator-Lr solution. Samples were barcoded using the Cell-ID 20-Plex Pd Barcoding Kit and combined with the Maxpar cell acquisition solution and EQ four-element calibration beads followed by acquisition on the Helios CyTOF system. Data were collected as flow cytometry standard (FCS) files, debarcoded and normalized to the acquired calibration bead signal.

Processing of Flow Cytometry Standard (FCS) files
The workflow of further processing of the FCS files is shown in Supplementary Fig. 1B and described in detail in Supplementary Materials and Methods. Briefly, normalized FCS files were analyzed by the Maxpar Pathsetter software (version 2.0; Fluidigm) for cleanup to select live cell events. To assess batch variability, unstimulated and cytokine-activated Veri-Cell reference standards (see Supplementary Materials and Methods) were compared between the batches. As expected [27], PMA and ionomycin stimulation resulted in consistent increases in pp38, pERK, pCREB, Ki-67, and pS6 and reductions in total IkBα ( Supplementary Fig. 1C). The markers SOX2, CD27, and CD147 were subsequently removed from the downstream analysis because of batch variations among patient samples. Additional processing steps including cell modeling (Supplementary Figure  1D) including criteria for inclusion of MM cells ( Supplementary Fig. 1E) and further data processing are described in Supplementary Materials and Methods.

Comparison to clinical metrics
The log2 fold change of marker expression was calculated to identify differences in the average marker abundance, proportion of cells with marker changes from clusters, and proportion of cells assigned to metaclusters for the 90 day treatment response (poor vs. good), sex (male vs. female), age ( > 60 vs. ≤ 60 years), type (relapsed refractory MM (RRMM) vs. newly diagnosed (NDMM), primary (hyperdiploidy vs. 11;14), TP53 deletion (present vs. absent), 1q gain (present vs. absent), MYC disruption (present vs. absent), monosomy 13 (present vs. absent), deletion 13q (present vs. absent), mSMART (high vs. standard), and S-phase (high vs. low). Statistical significance was determined using an unpaired t-test (ttest of independent means) between values for each comparison. The correlation coefficients between a pair of binary variables were calculated with Cramer's V method with signs (+/−) introduced from a Pearson's correlation test, which was used to estimate all of the remaining correlations. Violin plots were generated using the GraphPad Prism software (San Diego, CA).

Survival analysis
Overall survival (OS) was defined as the time from diagnosis to death from any cause or the last follow-up, with those alive censored at the date of the last follow-up. The APEX phase-3 clinical trial was used as a validation cohort. Pretreatment gene expression and clinical outcomes data on the APEX phase-3 clinical trial (039; n = 156) on relapsed and/or refractory myeloma patients were downloaded from the Gene expression omnibus (GEO) database (GSE9782) as Affymetrix HG-U133A/B gene probe set analysis dataset. The APEX data set (GSE9782; n = 264) contains gene expression data on phase-2 and phase-3 relapsed and/or refractory myeloma clinical trials, including APEX phase-3 trial (039), a APEX companion study (040), SUMMIT (025) and CREST phase-2 trials (024) [34]. Median gene intensity scores were calculated for the following genes that served as meta-cluster 1 gene signature: BCL-2, IKZF3, MYC, NFKBIA (IkBα), and PTPRC (CD45). A combined meta-cluster activity score was computed for each sample as the summation of the ascending ranks of expression of median gene intensities. Survival curves were estimated using Kaplan-Meier analysis and compared using the log-rank test. Cox regression analysis was performed to estimate the survival risk of the different phenotypic meta-clusters for overall survival. Statistical analyses were performed using BlueSky (Chicago, IL, USA), with a value of significance defined as P < 0.05.
CyTOF analysis recapitulates expected protein expression patterns in MM First, we analyzed our data and compared our findings to previously published data in MM. We compared the average normalized signal of each marker within each patient's bulk PC population across the dataset. Among the cell surface proteins, CD38, CD49d, CD138, and CD71 were the most abundant, while CD20, a late-stage B cell marker, was the least abundant, as previously reported [35] (Fig. 1A). Expression of CD56, a transmembrane glycoprotein known to be aberrantly expressed in MM [36], was present in 24 (49.0%) and absent in 25 (51.0%) samples (Fig. 1A, B). Consistent with previous reports [37,38], the average surface expression of the PC marker CD138 was lower than that of CD38; cryopreservation and sample processing have been reported to reduce CD138 surface expression. In addition, detection of surface CD38 protein is influenced by treatment with the CD38-targeting antibody daratumumab [39], interfering with CD38 identification 4-6 months after the last exposure [40][41][42]. Accordingly, the fourteen samples (28.6%) from patients who had   Fig. 1C, Supplementary Table 4).
In addition, we identified numerous proteins with higher abundance in RRMM compared to NDMM samples (Fig. 1B). The expression levels of BCL-2 (P < 0.01), clCasp3 (P = 0.04), Ki-67 (P = 0.04) and CD81 which has been associated with chemoresistant MM tumor cells [47] (P < 0.01), were significantly higher in RRMM than in NDMM samples (Fig. 1B) consistent with a more aggressive disease. These data using bulk protein expression demonstrate that CyTOF profiling recapitulates key protein features of MM as previously reported.
Thirteen commonly expressed phenotypic meta-clusters can be detected among patients Next, we determined whether our integrated cell surface and intracellular protein dataset could be used to detect clinically meaningful patterns of protein expression among subpopulations of MM tumors. First, to identify common phenotypic meta-clusters among our patient cohort, we used an unsupervised clustering analysis to define phenotypic clusters in each patient sample individually ( Fig. 2A and Supplementary Materials and Methods). This analysis revealed 11-29 phenotypic clusters per patient. A representative sample of RRMM with t(11;14) and TP53 deletion is shown in Fig. 2A.
To define the phenotypic clusters further, we next performed a differential abundance analysis to identify significant markers that differed in their expression between clusters in each patient (Fig. 2B). This analysis generated Z-scores for each marker showing the degree of directionality of the change in expression (increased expression, positive Z-scores and decreased in expression, negative Z-scores). Next, to determine if phenotypic clusters were shared across patients, individual patient clusters were divided into metaclusters using the Z-score data as inputs for unsupervised clustering to identify cluster-matched equivalents among all patients (Fig. 2C). This allowed us to characterize similar phenotypic clusters across all patients revealing 13 phenotypic meta-clusters among the 49 MM samples (Fig. 2D). Each phenotypic meta-cluster was defined by the most highly differentially expressed protein markers (absolute Z-score of > 5 and significant changes in > 50% of cells within the meta-cluster).
Phenotypic meta-cluster composition is correlated with clinical behavior Next, we asked whether the proportion of individual phenotypic meta-clusters within each patient correlated with clinical variables including tumor genetics, response to therapy and disease outcome (Fig. 3A). The fold change in the proportion of cells within a sample belonging to each meta-cluster is indicated in the heatmap and correlated to clinical variables (Fig. 3B). Tumors associated with a poor depth of response had a significantly reduced abundance of meta-cluster 1, which was characterized by elevated CD45 and low survival protein BCL-2 along with elevated IKZF3, MYC and IkBα, suggesting that phenotypic meta-cluster 1 may be responsive to MM therapy. Meta-cluster 13, characterized by increased expression of the proliferative marker Ki-67 and CD117 as well as the IKZF1 transcription factor and pRB (Fig. 2D,   Fig. 3 Comparison of the 13 phenotypic meta-clusters with clinical metrics, patient demographics, tumor genetics and disease outcome. A Overall schematic of workflow from subpopulation processing to clinical metric correlation. Meta-cluster 1 (red circle) was the only subpopulation associated with a significant difference in response to therapy and disease outcome. B Heatmap showing the log2 fold change in the proportion of cells within a sample belonging to each meta-cluster for different clinical and cohort metrics. Significant (P < 0.05) changes are highlighted in a black box. Cohort metrics defining response to treatment (poor including PD vs. good including PR and VGPR), sex (male vs. female), age (> 60 vs. ≤ 60 years), type (RRMM vs. NDMM)), primary genetic abnormality (hyperdiploidy vs. 11;14), TP53 deletion (present vs. absent), MYC disruption (present vs. absent), monosomy 13 (present vs. absent), deletion 13q (present vs. absent), mSMART score (high vs. standard), and S-phase value (high vs. low). C A comparison of OS (years) in patients with absence of meta-cluster 1 (red line) and in patients with presence of meta-cluster 1 (blue line). OS time (median) was 2.2 [95%CI: 1.6-3.9] years (n = 7) and 9.4 [95%CI: 3.6-not calculated] years (N = 17) in the 2 meta-clusters, respectively. D Distribution of high and standard mSMART (TP53 deletion and/or 1q gain) among the 7 patients without meta-cluster 1 and 17 patients with meta-cluster 1. Table 5) was significantly more prevalent in RRMM samples and in those with 1q gain, elevated mSMART and high S-phase (Fig. 3B). Meta-cluster 13 may be associated with disease relapse and suggests that patients with elevated meta-cluster 13 may benefit from anti-proliferative therapy. Further, samples associated with high mSMART also had a high proportion of meta-cluster 9, characterized by elevated cell surface transferrin receptor CD71, and pERK suggesting targeting CD71 may be a promising therapeutic strategy in these patients. Samples with high S-phase had an increased abundance of meta-cluster 2, characterized by elevated cell surface CD81 expression, another potential therapeutic target. Meta-cluster 8, characterized by the highest relative expression of CD38, was significantly reduced in RRMM samples and in those with TP53 deletion and high S-phase, possibly driven by daratumumab treatment.

Supplementary
Since CD38 was significantly lower in daratumumab-treated vs. daratumumab-naïve samples (Fig. 1C), we investigated whether differences in the phenotypic profiles of CD38-driven metaclusters were skewed due to daratumumab. We repeated clustering analysis and meta-cluster classification with the exclusion of CD38. As before, we determined significant markers with relative increases or decreases in expression in each cluster and compared the meta-clusters that had been analyzed with CD38 inclusion to those with CD38 exclusion (Supplementary Fig.  4A). In comparison to the previous analysis, in which 13 different phenotypic meta-clusters were identified, reanalysis without CD38 revealed 12 phenotypic meta-clusters. The expression characteristics of these 12 phenotypic meta-clusters were remarkably similar, suggesting that within these 12 meta-clusters, the influence of CD38 on meta-cluster assignment was minimal. Meta-cluster 8, characterized in the original analysis as having the highest CD38 expression (Supplementary Table 5), was the only meta-cluster lost in the reanalysis suggesting that CD38 played a significant role in the classification of meta-cluster 8. Approximately 11% of the cells belonging to meta-cluster 8 in the previous analysis were redistributed to other meta-clusters ( Supplementary Fig. 4A). Similar to the original analysis, metacluster 13, with increased Ki-67 and CD117, remained significantly more prevalent in samples with 1q gain and elevated mSMART. Samples associated with a poor depth of response retained their reduced abundance of meta-cluster 1 in the reanalysis (Supplementary Fig. 4B).
We next compared whether the presence or absence of each of the 13 phenotypic meta-clusters within 6 months of MM diagnosis was predictive of overall survival (OS) in our patient cohort. The only phenotypic meta-cluster associated with a significant difference in OS was meta-cluster 1 (Fig. 3C). Among patients with an absence of phenotypic meta-cluster 1, OS was decreased compared to patients with the presence of phenotypic meta-cluster 1 (2.2 years vs. 9.4 years, P = 0.017) characterized by increased CD45 and reduced BCL-2 expression. Consistent with these findings, the presence of phenotypic meta-cluster 1 was associated with improved OS (RR 0.21 (0.05-0.84), P = 0.027) (Supplementary Table  6). A similar distribution of mSMART risk status was observed in association with the presence (52.9%) or absence (57.1%) of phenotypic meta-cluster 1 (P = 0.67) (Fig. 3D). In a multivariate model including the presence of meta-cluster 1 and standard risk mSMART status, meta-cluster 1 retained its trend toward improved OS (RR 0.28 (0.07-1.08), P = 0.063) (Supplementary Table 6).
We sought to validate the association of meta-cluster 1 with a survival benefit in a larger patient cohort. Since an external CyTOF dataset is not available, we utilized gene expression and clinical outcome data from the APEX phase-3 MM clinical trial. Median gene intensity scores were calculated for the genes that defined the meta-cluster 1 protein signature: high PTPRC (CD45) and low BCL-2 as well as high NFKBIA (IkBα), MYC and and IKZF3. A combined meta-cluster 1 activity score was computed for each sample as the summation of the ascending ranks of expression of median gene intensities. To determine whether the meta-cluster 1 activity score was associated with clinical outcomes, we performed survival analysis between patients with the highest vs. the lowest meta-cluster 1 score. Patients with a low metacluster 1 score were associated with poorer OS in comparison to patients with a high meta-cluster 1 score (Log-rank hazards ratio = 2.052 (95% CI 1.121-3.754; p = 0.018) (Fig. 4). Further, the odds of no response (NR) vs. response (R) between the patients with the lowest meta-cluster 1 score vs. the highest meta-cluster 1 score was 2.6786 (PGx responder status vs. meta-cluster 1 score; Odds ratio (OR) = 2.6786; P = 0.0338; 95% CI 1.0781 to 6.6549), indicating a low meta-cluster 1 score was associated with > 2.5fold higher cancer progression in MM patients compared to a high meta-cluster 1 score.

DISCUSSION
Here we used single-cell CyTOF to characterize MM tumor heterogeneity among 49 patient samples. Our goal was to detect shared MM subpopulations associated with treatment response and disease outcome across all patient samples. We employed differential expression analysis to group similar clusters within each patient sample so that we could compare phenotypic metaclusters among all patients. Using a random forest model, we identified cells with similar marker expression patterns, allowing us to identify 13 unique phenotypic meta-clusters. To the best of our knowledge, this is the first study to evaluate subclonal heterogeneity, characterized by the analysis of~40 phenotypic cell surface and intracellular markers using primary MM samples. Despite significant heterogeneity, we found common patterns of subclonal protein profiles that were shared across patients and are correlated with clinical behavior. Evaluation of tumor proteomic profiles may provide an opportunity for the identification of novel protein-based biomarkers associated with disease state and Of the 13 phenotypic meta-clusters identified, meta-cluster 1 was the only meta-cluster that had a significant impact on response to therapy and overall disease outcomes. This impact was validated using the APEX 3 clinical trial. Meta-cluster 1 was primarily characterized by a subpopulation of MM cells with increased CD45 and reduced BCL-2 expression. Reduced abundance of meta-cluster 1 was associated with poor treatment response and survival. This observation is consistent with previous studies demonstrating that patients with CD45 negative MM had poorer OS compared to patients with CD45 positive MM, suggesting that the CD45 negative MM phenotype may be a marker of progressive disease [48,49]. CD45 low or negative MM cells have also been identified to have a higher engraftment ability in mice compared to CD45 positive MM cells [50]. Of interest, the CD45 positive MM cells have been previously associated with low expression of the pro-survival protein BCL-2 along with an increased susceptibility to apoptosis [51], suggesting that phenotypic meta-cluster 1 may be more sensitive to treatment-related cell death. Our work also revealed increased BCL-2 expression in RRMM, a finding that has previously not been reported and may suggest an increased sensitivity to venetoclax in this patient cohort. Further, since we show that BCL-2 expression correlates with t(11;14), the increased BCL-2 expression in RRMM was not a reflection of a higher frequency of t (11;14) in the RRMM group (NDMM: 35% t(11;14) vs. RRMM: 28% t(11;14)).
Consistent with previous studies [41,52], we found that the monoclonal CD38 antibody daratumumab interfered with the CD38 antibody used in our CyTOF studies, precluding the reliable use of CD38 for initial gating approaches and in the downstream analysis of patients treated with daratumumab within 6 months prior to sample collection. Remarkably, the exclusion of CD38 in the reanalysis of our data appeared to primarily impact the identification of phenotypic meta-cluster 8, characterized by the highest CD38 expression, and the remaining 12 clusters maintained their phenotype identity independent of CD38 expression. Future studies using alternative CD38 epitopes that do not compete with daratumumab should be conducted. Although the overall detection of surface CD38 was lower in samples obtained from patients treated with daratumumab than in untreated patients (Fig. 1C), two patients represented by the highest two data points in the treatment cohort did not have reduced CD38 detection 2 and 4-months post treatment. Whether this was a result of altered CD38 expression [40,41] or a lack of antibody interference by daratumumab remains unknown [41].
This study is limited by its retrospective nature owing to the use of cryopreserved samples for CyTOF analysis and heterogeneity in treatment regimens in the RRMM group (Supplementary Table 4). Despite these limitations, we showed that MM samples display a significant degree of complexity in protein expression, but that subclonal structure is shared among patients. We identified 13 unique phenotypic meta-clusters across the 49 MM samples. Several of the meta-clusters were associated with unique clinical features and disease subtypes. The marker proteins of each phenotypic meta-cluster may serve as therapeutic targets or biomarkers which could be evaluated in future studies.

DATA AVAILABILITY
Requests regarding data availability not already shared in the supplementary material should be made to the corresponding author.